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We study the microscopic structure and the stationary propagation velocity of (l + l)-dimensional 
solid-on-solid interfaces in an Ising lattice-gas model, which are driven far from equilibrium by an 
applied force, such as a magnetic field or a difference in (electro)chemical potential. We use an 
analytic nonlinear-response approximation [P. A. Rikvold and M. Kolesik, J. Stat. Phys. 100, 377 
(2000)] together with kinetic Monte Carlo simulations. Here we consider interfaces that move under 
Arrhenius dynamics, which include a microscopic energy barrier between the allowed Ising/lattice- 
gas states. Two different dynamics are studied: the standard one-step dynamics (OSD) [H. C. 
Kang and W. Weinberg, J. Chem. Phys. 90, 2824 (1992)] and the two-step transition-dynamics 
approximation (TDA) [T. Ala-Nissila, J. Kjoll, and S. C. Ying, Phys. Rev. B 46, 846 (1992)]. In 
the OSD the effects of the applied force and the interaction energies in the model factorize in the 
transition rates (soft dynamics), while in the TDA such factorization is not possible (hard dynamics). 
In full agreement with previous general theoretical results we find that the local interface width 
under the TDA increases dramatically with the applied force. In contrast, the interface structure 
with the OSD is only weakly influenced by the force, in qualitative agreement with the theoretical 
expectations. Results are also obtained for the force-dependence and anisotropy of the interface 
velocity, which also show differences in good agreement with the theoretical expectations for the 
differences between soft and hard dynamics. Our results confirm that different stochastic interface 
dynamics that all obey detailed balance and the same conservation laws nevertheless can lead to 
radically different interface responses to an applied force. 

PACS numbers: 68.35.Ct 75.60. Jk 68.43.Hn 05.10.Ln 



I. INTRODUCTION 



To understand the properties of materials one has to understand surfaces and interfaces: materials interact with 
their environment through their surfaces, and material properties are profoundly influenced by internal interfaces. 
Processes related to surfaces and interfaces therefore play a critical role in nature and in a variety of technological 
applications, such as electronic, magnetic, and optical devices, sensors, catalysts, coatings, and many other industrial 
systems and processes. 0, Q Surfaces also play a vital role in biology and medicine, since most biological reactions 
occur at surfaces and interfaces. The applications of surface science in medicine range from the growth of bio- 
compatible surfaces for tissue cultures, through medical implants, to the design of innumerable medical devices. 0,3 
It is therefore crucial to understand the fundamental processes occurring at surfaces and interfaces in order to provide 
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means of controlling and manipulating such systems. The large-scale properties of growing interfaces have been the 
object of an enormous amount of work in recent years, 0] but much less attention has been paid to interfacial 
structure on a microscopic scale. This is unfortunate because many important interface properties, such as mobility 
and catalytic and chemical activity, are largely determined by the microscopic interface structure. The microscopic 
structure limits the interfacial propagation velocity under an external driving force, such as an external field for a 
magnetic or dielectric domain wall, or a difference in chemical potential between the bulk phases for a crystal surface. 

Since the detailed microscopic mechanism of the interface motion is often not known, one standard way to gain 
insight in the process is by constructing a stochastic model that mimics its essential features. Dynamics that conserve 
the order parameter, or ones that do not, must be chosen according to the physical characteristics of the system. Once 
this is decided, there are several dynamics in each category to choose among. It is well established that structures 
arising from different dynamics that obey detailed balance and respect the same conservation laws exhibit universal 
asymptotic large-scale features. However, recent studies 0,H> show that there are important differences between 

the micro structures of field-driven interfaces obtained with different dynamics, and that these differences very signifi- 
cantly influence important interface properties such as mobility. A mean-field, nonlinear-response theory developed in 
Refs. 0,0 indicates that there are two different classes of nonconservative dynamics that lead to significantly different 
surface microstructures. These are the soft dynamics in which the single-site transition rates can be factorized into 
one part that depends only on the applied force and a second part that depends only on the interaction energies, 
and the hard dynamics for which this factorization is not possible. jTH | By this classification, the widely used Glauber 
and Metropolis dynamics are hard. Soft dynamics are appropriate for solidification or adsorption problems where the 
driving force is a chemical-potential difference. [T^. Il3| 

In previous papers^, @ two of us introduced a dynamic mean-field approximation for the microstructure, based 
on the Burton-Cabrera- Frank solid-on-solid (SOS) approximation. For soft dynamics, interface structures should 
remain independent of the applied force, while there should be a clear dependence on the force for hard dynamics. 
Monte Carlo (MC) simulation results with the soft GlauberQ and the standard hard Glauber dynamics[a |j| 0] 
confirm these predictions. Moreover, SOS surfaces generated with the hard Glauber dynamics exhibit a skewness that 
is absent for the soft Glauber dynamics. Since the mean- field results depend on the absence of short-range correlations 
along the interfaces, it remains an open question to which extent these characteristics are shared by all soft and hard 
dynamics. 

For the present study we have considered two dynamics that include a local energy barrier representing a transition 
state inserted between individual Ising or lattice-gas states. Such Arrhenius dynamics, as they are often called, are 
appropriate in kinetic MC simulations of discrete Ising or lattice-gas models in which the discrete states serve as 
approximations for high-probability configurations in an underlying continuous potential. [TH 1 1 6l| Examples are the 
study of diffusion in a lattice-gas model in a continuous corrugation potential. |l7l| the relaxation process from the 
high-spin state of molecular bistable solid s. Il8l ITfll l2(ij | or the approximation of a continuous spin model with strong 
uniaxial anisotropy by an Ising model. [2ll l22l | When applied to kinetic Ising lattice-gas models, Arrhenius dynamics 
give nucleation rates quite different from the ones given by the standard Glauber dynamics; another warning about 
the importance of choosing the right dynamics for specific physical or chemical systems. 

The two Arrhenius dynamics that we consider here are the common one-step-dynamics fOSD)[24Ll25j] and the two- 
step transition dynamic approximation fTDA).[l^.l26l| The OSD dynamics are soft and the TDA dynamics are hard, 
according to the definitions given above. Their transition rates are defined in Sec. II. The OSD dynamics are commonly 
used in studies of adsorption, such as the electrosorption of halides on single-crystal silver electrodes. (23, 123, 123 
Among the experimental systems that have been described with the TDA dynamics are the diffusion of H atoms on 
single-crystal tungsten surfaces. [26j 

SOS surfaces belong to the Kardar-Parisi-Zhang (KPZ) dynamic universality class . [BL l30| in which the macroscopic, 
stationary distribution for flat, moving interfaces is Gaussian, corresponding to a random walk with independent 
increments. Nevertheless, the step heights in several discrete models in this class are correlated at short distances. |3ll 
[3^ In the mean-field approximation used here, these short-range correlations are ignored. The resulting discrepancies 
will be apparent when we compare the theoretical results with kinetic MC simulations. 

The remainder of this paper is organized as follows. In Sec. ^ we introduce the SOS interface model and give the 
transition rates for the TDA and OSD dynamics. In Sec. lIIII we summarize the mean-field approximation for the time 
evolution of the single-step probability density function (pdf), as well as its stationary form. We also give expressions 
for the class populations and interface velocity in terms of the applied force, the temperature, and the angle of the 
interface relative to the lattice axes. In Sec. IIVI we compare simulations and analytical predictions for the detailed 
stationary interface structure, including the asymmetry of the simulated nonequilibrium interfaces. A summary and 
conclusion are provided in Sec. 
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II. MODEL AND DYNAMICS 



The SOS interfaces are described by the nearest-neighbor S = 1/2 Ising Hamiltonian with anisotropic, ferromagnetic 
interactions J x and J y in the x and y direction, respectively: 

T~i ^ s x,y {JxS x +l,y + JySx,y+l + H), (1) 

x,y 

where s XtV = ±1, J2 X y runs over ai l sites, and the applied field H is the driving force. The interface is introduced 
by fixing s XjV — +1 and —1 for large negative and positive y, respectively. Without loss of generality we take H > 0, 
such that the interface on average moves in the positive y direction. This Ising model is equivalent to a lattice-gas 
model with local occupation variables c x , y £ {0, 1}.[33|,|34[ Specifically, we identify s = +1 with c = 1 (occupied or 
"solid") and s = — 1 with c = (empty or "fluid"). The interactions in the Ising model, J a , are related to the ones 
in the lattice-gas model, (j> a , as J a — </> Q /4, and the applied field is related to the lattice-gas chemical potential fi 
as H = (/i — fj,o)/2, where /io = —^{Jx + Jy) = — {4>x + 4>y) is the coexistence value of //. Here we will use Ising or 
lattice- gas language interchangeably as we feel it makes a particular aspect of the discussion clearer. 

The SOS model considers an interface in a lattice-gas or S = 1/2 Ising system on a square lattice of unit lattice 
constant as a single- valued integer function h(x) of the x-coordinate, with steps 5(x) = h{x + 1/2) — h(x — 1/2) at 
integer values of x. A typical SOS interface configuration is shown in Fig. ^ I n this paper the two possible states 
of the site (x,y) are denoted by the two "spin" values s x . y — ±1. (In order that the step positions and the interface 
heights be integer as stated above, we place the spins at odd half-integer values of x and y, i.e., at the centers of the 
unit cells separated by dotted lines in Fig.^) 

The interface will be made to evolve under two different dynamics that contain a microscopic energy barrier 
against individual spin flips. These are the transition dynamics approximation fTDA) pH l26| and the so-called one- 
step dynamics (OSD). 0,IH| The barrier represents a transition state which is inserted between the states allowed 
in the Hamiltonian, such as a saddle point in a corrugation potential for particle diffusion. pH, fill |2(| or a high 
energy associated with a transitional spin state that is not along one of the two directions allowed by the Ising 
Hamiltonian. plli^ 

Here we express the transition-state energy by the following approximation. [Hi Izil I2H |26 | 



E T = ^L+U, (2) 

where Ei and Ef are the initial and final energies, and U is the bare, microscopic energy barrier. In electro- 
chemical applications, such as electron- or ion-transfer reactions, this corresponds to the symmetric Butler- Volmer 
approximation. 35j The construction corresponding to Eq. @ is illustrated in Fig. [5] 

Both the TDA and the OSD are single-spin-flip (nonconservative) dynamics that satisfy detailed balance. This 
ensures the approach to equilibrium, which in this case is a uniformly positive phase with the interface pushed off 
to positive infinity. Such dynamics are defined by a single-spin transition rate, W(s x . y — > — s x . y ) = W((3AE, (3U). 
Here (3 is the inverse of the temperature T (Boltzmann's constant is taken as unity), AE — Ef — Ei is the en- 
ergy change corresponding to a successful spin flip, and U determines the energy barrier between the two states 
through Eq. J2J). The detailed-balance condition (valid for transitions between allowed states) is expressed as 
W(f3AE, (3U)/W(-f3AE, /3U) = e^ 0AE , where the right-hand side is independent of U. 

The transition rates for the TDA and the OSD dynamics with the transition-state energy Et given by Eq. © are 
defined as [3 EHH 

Wtda = (3) 

TDA 1 + exp[/?(£ T - Ei)] 1 + exp\p(E f - E T )} { ' 

and 0,0 

Wosu = eM~P(E T ~ Ei)] = exp(-pU) exp[-/3A£/2], (4) 

respectively. 

Notice that the TDA transition rate cannot be factorized into one part that depends only on the interaction energy 
and another that depends only on the applied field; thus it belongs to the class of dynamics defined as hard.0] The 
OSD dynamics can be factorized this way and thus is classified as soft.0] Another important difference between the 
TDA and the OSD is that Wtda is restricted to be between and 1, while there is no upper bound on Wosd- As a 
consequence, the same difference is observed for the propagation velocities. One may however question the physical 
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realism of cases in which the transition-state energy Et is below the initial energy Ei so that Wosd > 1. In order to 
preserve the SOS configuration at all times, flips are allowed only at sites which have exactly one broken bond in the 
y direction. 

With the Ising Hamiltonian there are only a finite number of different values of AE. The spins can therefore be 
divided into classes. |36L 1371 1381 139] labeled by the spin value s and the number of broken bonds between the spin and 
its nearest neighbors in the x and y directions, j and k, respectively. The ten spin classes consistent with the SOS 
model are denoted jks with j £ {0, 1, 2} and k £ {0, 1}. They are shown in Fig. H an d listed in Tablc[lJ 

In the SOS model the heights of the individual steps are assumed to be statistically independent and identically 
distributed. The step-height probability density function (pdf) is given by the interaction energy corresponding to 
the \S(x)\ broken J^-bonds between spins in the columns centered at (x — 1/2) and (x + 1/2) as 



p[6{x)} = Z{<j))- l X Wx) \ e^ 6 ^ 



(5) 



The factor X determines the width of the pdf, and "f(4>) is a Lagrange multiplier which maintains the mean step 
height at an ^-independent value, (S(x)) — tan0, where <fi is the overall angle between the interface and the x axis. 



In equilibrium, X is simply the Boltzmann factor, e 



-2/3J x 



The partition function is, 



z{4>) = 



£ 

S=— oo 



l-X 2 
1 - 2X cosh 7 (0) + X 2 



(6) 



where 7(</>) is given by 



1 + X 2 )tan< 



'1-X 2 rtan 2 



AX< 



1/2 



(see details in Refs. 



2X{1 + tan</>) 

Simple results are obtained for <j) = 0, which yields 7(0) = and 
Z(0) = (l+X)/(1-X), 



(7) 



(8) 



and for = 45° (see Ref. i). 

The mean spin-class populations, (n(j'fcs)), are all obtained from the product of the independent pdfs for S(x) and 
S(x+1). Symmetry of p[S(x)] under the transformation (x,4>,5) — > (— x, — cf), — S) ensures that (n(jk—)) = (n(jk+)} 
for all j and k. Numerical results illustrating the breakdown of this up/down symmetry for large H are discussed 
in Sec. IIVI The general expressions for the class populations are given in the first column of Table ITU details of the 
calculation can be found in Ref. |jj The results for each of the dynamics are obtained by substituting their respective 
values of X , which will be calculated in the next Section. 

Whenever a spin at the interface flips from —1 to +1, the corresponding column of the interface advances by one 
lattice constant in the y direction. Conversely, a column (not necessarily the same one) recedes by one lattice constant 



when a spin at the interface flips from +1 to —1. The mean velocity of the interface in the y-direction, 



is the 



difference between the rates of forward and backward steps, averaged over the whole interface. The energy changes 
corresponding to the flips are given in the third column in Table [TJ The sum over x that gives the average velocity can 
be rearranged into a sum over the spin classes of transition rates, weighted by the average class populations. Since 
the spin-class populations on both sides of the interface are equal in this approximation, the contribution to (v y ) from 



sites in the classes jk— and jk+ therefore becomes 

(v y (jk)) = W((3AE(jk-),pU)-W(f3AE(jk+),f3U) . 
The mean propagation velocity perpendicular to the interface becomes 

{v±(T,H,(l>)) = cos(4>){v y ) = cos(<t>)J2Wjks)}{v y (jk)) 



(9) 



(10) 



where the sum runs over the classes included in Table [TJ] It has been shown in Ref. that Eq. (|1L)|) reduces to the 
results for the single-step [3^. liol l4ll |42| and the polynuclear growth [iol liM lil] models at low temperatures for large 
and small 0, respectively. 
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III. NONLINEAR RESPONSE 



With X = e~ 2/3Jx , i.e, independent of H, the results in Table HH correspond to a linear-response approximation 
for the velocity. In previous paperspi Q an expression for a field-dependent X(T,H) was obtained, based on a 
dynamic mean-field approximation for the equation of motion for the single-step pdf together with a detailed-balance 
argument for the stationary state. This improved non-linear response approximation gives (see Ref. for details of 
the calculation), 

X(T m = -2 PJx { e-^ H W[f3{-2H 4.U (3U] + e^ H W{f3(2H 4J a ), f3U] \ 1/2 
1 ' ] I W[f3(-2H-4J x ),/3U]+W[P(2H-U x ),pU] J ' 1 ' 

which is independent of 7(0). Here W(j3AE, f3U) are the transition rates associated with the reversal of a single spin. 
The values of AE, Et — Ei, and Ef — Et are given in Table [I] for the different spin classes. 

Equation shows that X(T,H) depends on the specific dynamics, except for H = 0, where it reduces to its 
equilibrium value, X(T,0) — e~ 2/3Jx . It is easy to see from the equation that for soft dynamics, where the field and 
interaction contributions to the transition rates factorize, the ii-dependence in Eq. (|llf> cancels out. In Ref. |JJ it 
was demonstrated that the soft Glauber dynamics yields an SOS interface that is identical to the equilibrium SOS 
interface at H = and the same temperature, regardless of the value of H. Hard dynamics, such as the standard 
Glauber and Metropolis dynamics and the TDA, lead to a nontrivial field dependence in X . 

Inserting the transition rates corresponding to the TDA and the OSD dynamics, defined by Eq. © and Eq. 
respectively, into Eq. pi[). we explicitly get 



x (T m _ -20 Jx j e 2PJ * cosh(2/3ff) + er 2 ^ + 2coshQ3£/) cosh(/3ff) \ 1/2 
TDM ' 1 \ e- 2 ? J * cosh(2/3£0 + e 2 ^ + 2cosh(/3C/) cosh(/3iT) J 1 ' 

and 

^osd(T, H) = e- 2 P J * = X(T, 0). (13) 

We note that Atda is similar, but not identical, to the one for the standard Glauber dynamics, Eq. (18) of Ref. 0. 
The spin-class populations listed in Table Irfl can now be calculated explicitly for each of the dynamics by replacing X 
with its corresponding value. The expressions for the contributions to the mean velocity in the y direction, Eq. J^J, 
for each class in the TDA and the OSD dynamics are given in the third and fourth columns of Table [H] respectively. 

In the next Section we show that the nonlinear-response approximation gives good agreement with MC simulations 
of driven, flat SOS interfaces evolving under the TDA and OSD dynamics for a wide range of fields and temperatures. 



IV. COMPARISON WITH MONTE CARLO SIMULATIONS 



We calculated the step-height distributions, propagation velocities, and spin-class populations, analytically and by 
kinetic MC simulations, for both the TDA and the OSD dynamics in the isotro pic case, J x = J y = J. The details 
of our particular implementation of the n-fold way rejection-free MC algorithm|37l l3gj are essentially the same as 
described in Ref. 0, except for two points. The first is that only transitions from the classes with one broken y-bond 
(k = 1) are allowed, so as to preserve the SOS interface structure. The second difference is that the present code 
uses continuous time^j^ to accommodate the large transition rates that are possible with the OSD dynamics. By 
keeping only the interface sites in memory, the algorithm is not subject to any size restriction in the y direction, and 
simulations can be carried out for arbitrarily long times. 

The numerical results presented here are based on MC simulations mostly at the two temperatures, T = 0.2T C 
and 0.6T C (T c = —2 J/ ln(V2— 1) ~ 2.269 J is the critical temperature for the isotropic, square-lattice Ising modelpq'). 
with L x = 10 000 and fixed 4> between and 45°. The microscopic transition barrier, U (see Eq.(2)), is chosen to be 
0.5 J. This is the same value used in a previous study of nucleation with the OSD and TDA dynamics. [23 From Eq. (@J 
it is clear that for the OSD dynamics U only appears in a temperature-dependent scaling factor in the transition rate. 
It thus has no influence on the interface structure. For the TDA dynamics, on the other hand, an increase in U leads 
to a decrease in the local interface width and consequently in the propagation velocity. Also, the observed skewness 
increases, suggesting increasing short-range correlations between the step heights. These results are discussed in detail 
in a separate paper. [43 

In order to ensure stationarity we ran the simulation for 50 000 n-fold way updates per updatable spin (UPS) before 
taking any measurements. Stationary class populations and interface velocities were averaged over 50 000 UPS. For 
the stronger fields at T = 0.2T C we used ten times as many UPS. Adequate statistics for one- and two-step pdfs were 
ensured by the large L x . 
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A. Stationary single-step probability densities 



Stationary single-step pdfs were obtained by MC simulation at T — 0.2T C and 0.6T C for (ft — and several values 
of H. The simulation data and the theoretical results for p[S] are shown in Fig. Eland Fig. 01 for the TDA and the 
OSD dynamics, respectively. The theoretical results are calculated with Eq. JSJ, with X(T,H) from Eq. (|12(l for the 
TDA dynamics and from Eq. (|13(l for the OSD dynamics. For the TDA dynamics, Fig. |21 the agreement is excellent 
at the higher temperature for all the values of H analyzed (up to H/J — 4). However, at the lower temperature 
the agreement is not very good for fields above H/J = 2.5. For the OSD dynamics, Fig. 01 shows that, contrary to 
the theoretical mean-field results, p(S) depends somewhat on H. This dependence is stronger for small fields and 
at the lower temperature. However, although not absent as expected from the mean-field approximation, the field 
dependence is much weaker than for the TDA dynamics. 

Another way to compare the analytical and simulation results is by calculating (\5\) by summation of Eq. 
{\8\) = 2X/ (1 - X 2 ), with X from Eq. O for the TDA dynamics or X from Eq. (JTSJ (independent of H) for the 
OSD, and comparing these values with the simulated ones. The simulation values for (|<5|) can be obtained in two 
ways: directly by summation over the numerically obtained pdf, and also from the probability of zero step height as 
{\5\) = {p[0] _1 — p[0]} /2. This latter expression is obtained by observing that p[0] = (1 — X)/(l + X) for (ft — 0, then 
solving this for X, and inserting the result in the above equation for (|<5|) in terms of X. The results are shown in 
Fig-Elfor both dynamics, for (ft = at T = 0.2T C and 0.6T C calculated theoretically (solid lines) and by MC simulation 
(symbols). The agreement between the simulation and theoretical results for the TDA is excellent except at the low 
temperature, where for intermediate fields a slight deviation can be seen. The results are similar to those obtained 
with the standard Glauber dynamics. (See Fig. 5(a) of Ref. $,.) Again the results for the OSD dynamics show that, 
contrary to the theoretical prediction, there is a clear, albeit weak, dependence of the step height on the field. The 
theoretical and simulation results for the OSD dynamics only coincide at H/J = 0. 

The difference between the two dynamics is evident: the step heights for the OSD dynamics are weakly dependent 
on H, particularly for low values oi H, in contrast with the very strong H dependence obtained for the TDA dynamics. 
This behavior is typical for differences expected between soft and hard dynamics. 



B. Stationary interface velocities 

In this Section we compare the interface velocities obtained for the TDA and the OSD dynamics. In each case the 
velocities are calculated with the analytical approximation, Eq. I|1U[1 . and by simulations. Fig.||Ja) and Fig.||Jb) show 
the normal velocity vs H for (ft — for the TDA and the OSD dynamics, respectively. There is excellent agreement 
between the MC results and the nonlinear-response theory for the TDA dynamics, except for a slight disagreement 
seen at 0.2T C between H/J = 1.5 and 2.5. The results are very similar to those obtained with the standard Glauber 
dynamics. (See Fig. 6 of Ref. la) However, for the OSD dynamics at the lower temperature, the nonlinear-response 
approximation underestimates the velocity, especially at higher fields. One of the main differences between the two 
dynamics is clearly seen in Fig. the velocity is bounded by unity for the TDA while it increases exponentially with 
H for the OSD. 

The dependence of the normal velocity on the tilt angle (ft is shown in Fig. [7| and Fig. |S1 for the TDA and the 
OSD dynamics, respectively. We show results for several values of H/J at T = 0.2T C and T = 0.6T C . For the TDA 
dynamics, the agreement between the theoretical results and the simulations is excellent. The results are qualitatively 
similar to those obtained with the standard Glauber dynamics. (See Fig. 7 of Ref. 0.) For the OSD, the agreement 
between theory and simulation is excellent at the higher temperature. However, at the lower temperature and higher 
fields the agreement is only good at large values of (ft. 

In both cases it is seen that at T = 0.2T C in weak fields the velocity increases with (ft, in agreement with the 
polynuclear growth model at small angles and the single-step model for larger angles. However, for strong fields the 
TDA dynamics change gradually to the reverse anisotropy of Eden-type models. |H,E3 No such change is observed 
for the OSD dynamic. At T = 0.6T C , on the other hand, both models behave very similarly. The velocity is nearly 
isotropic for weaker fields, while becoming Eden-like for stronger fields. 

The temperature dependence of the normal interface velocity is shown in Fig. Elfor several values of H/J. For 
the TDA dynamics, Fig. |^a), the agreement between the simulations and the analytical results is excellent almost 
everywhere, except for a small discrepancy at intermediate T and H . This discrepancy is expected from the results 
shown in Fig. This figure shows that at low T, the velocity changes steeply from zero to unity at some value 
between H/J = 2 and H/J = 2.5, developing a step discontinuity in H at T = 0. The results are qualitatively similar 
to those obtained with the standard Glauber dynamics. (See Fig. 8 of Ref.0) For the OSD, Fig.^b), the theoretical 
and simulation results agree only when the temperature is higher than a minimum value that increases as the field 
increases. In this case the velocity also has a step discontinuity at T — 0: for H/J below 1.5 the velocity goes to zero. 
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For stronger fields the velocity at T = increases dramatically with H and decreases very rapidly as the temperature 
increases. As can be seen from Eq. JTHl and Table CU the mean-field theory predicts that the contributions to the 
propagation velocity from each of the classes at low T is proportional to exp [—/3 (2J X + U — H)]. Thus, there is a 
discontinuity at T = for H = 2J X + U, i.e., for H/J = 2.5 for our selection of U = 0.5 J. Beyond this value of H, 
the velocity diverges to infinity as T — > 0. Such a diver genc e at T = is not present for the soft Glauber dynamics, 
where the velocity vanishes at T — for all values of H .jlfij 



C. Spin-class populations and skewness 



Since the analytical predictions for the class populations are based on the assumption that different steps are 
statistically independent, a comparison with the simulation results gives a way of testing this assumption. The six 
mean class populations - (n(01s)), (n(lls)), and (n(21s)) with s = ±1 - for = at T = 0.2T C and 0.6T C are shown 
vs H in Fig. 1101 and Fig.^]f or the TDA and the OSD dynamics, respectively. For the TDA at both temperatures, 
Fig. 1101 the analytical approximations follow the average of the populations for s = +1 and s = — 1 quite well. 
However, at intermediate fields the simulations show that the population in front of the surface (s = —1) is quite 
different from the one behind it (s — +1). The mean-field approximation seems to reproduce better the population 
behind the surface. The results are qualitatively similar to those obtained with the standard Glauber dynamics. (See 
Fig. 9 of Ref.0.) For the OSD, Fig. ^2 the mean-field approximation predicts that the mean class populations should 
be independent of H (since X is independent of H) , while the simulation indicates a weak 7J-dependence. For small 
fields there is a clear dependence of the population on the field, but as the field increases, the populations tend to 
fixed values independent of the field. The iJ-dependence is consistent with the results for (\5\) and X(T,H), shown 
in Fig. EOT.) . 

The short-range correlations between neighboring steps are responsible for the skewness between the spin popula- 
tions on the leading and trailing edges of the interface that appears in the simulation results. This phenomenon is 
commonly observed in driven interfaces. It occurs even when the long-range correlations vanish as they do for inter- 
faces in the KPZ dynamic universality class, to which the present models belong for all finite, nonzero values of H. 
Skewness has also been observed in several other SOS-type models, such as the body-centered SOS model studied by 
Neergaard and den Niis.|3l| the model for step propagation on crystal surfaces with a kink-Ehrlich-Schwoebel barrier 
studied by Pierre-Louis et al..|50j and a model for the local time horizon in parallel kinetic MC simulations studied 
by Korniss et al.[32^ No skewness was observed for the SOS model with the soft Glauber dynamics. However, a 
small skewness was observed for the Ising model (whose interfaces include bubbles and overhangs) with soft Glauber 
dynamics (about two orders of magnitude smaller than the skewness observed for the hard Glauber dynamics). |l0j| 
The correlations associated with the skewness generally lead to a broadening of protrusions on the leading edge ( "hill- 
tops"), while those on the trailing edge ("valley bottoms") are sharpened. [31| or the other way around. |32j In terms 
of spin-class populations, the former corresponds to (n(21— )) > (n(21+)) and (n(ll+)) > (n(ll— )). The relative 
skewness can therefore be quantified by the two functions. |3l| 

_ _ <n(21-)) - <n(21+)) 



(n(21-)) + <n(21+)) 



and! 



(n(ll+))-(n(ll-)> 

(n(ll+)> + <n(ll-)) v ' 

These two skewness parameters are shown together in Fig. I12f at for the TDA, and in Fig. I12f bt for the OSD. For 
both dynamics the relative skewness is seen to be considerably stronger at the lower temperature. The temperature 
dependence is especially pronounced for p, due to the low concentration of sites in the class 21+ at low temperatures. 
The OSD results contrast with previous results that showed that an SOS interface with soft Glauber dynamics does 
not present skewness. |9j In the present case, the SOS interface evolving under the soft OSD dynamics presents stronger 
relative skewness at lower temperature than the one evolving under the hard TDA dynamics. 

Another way to visualize the skewness is to consider the joint two-step pdf, p [S(x), 5(x + 1)]. Logarithmic contour 
plots of this quantity for both dynamics at different values of H, for <f> = at T = 0.6T C , are shown in Fig. 1131 and 
Fig. [21 It is clearly seen that in both dynamics the contours change with H . For H=0 a symmetric diamond shape 
with equidistant contours indicates statistical independence with single-step pdfs given by Eq. (JSJ. This equilibrium 
result is correctly observed in both the TDA dynamics, Fig. I13f a~). and the OSD dynamics, Fig. I14f a1. However, for 
nonzero fields the dynamics show different behavior. For the TDA, at stronger fields, the shapes are convex in the 
second quadrant [6(x) < 0, S(x + 1) > 0] Fig. I13f ct and (d). This shape indicates that large negative S(x) tend 
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to be followed by large positive 5(x + 1) (sharp valleys). For weak fields, the TDA shows a weak skewness of the 
opposite sign, as seen in Fig. I13f b). For the OSD dynamics, for all nonzero fields, the shapes are concave in the 
second quadrant, Fig. H^T bl. (c), and (d). This shape indicates that large negative 5(x) tend to be followed by smaller 
positive S(x + 1) (rounded valleys). For both dynamics the field dependence in the fourth quadrant, corresponding to 
the shape of " hilltops," is much weaker that in the second quadrant. As expected, the contour plots for interfaces 
with <j) — are always symmetric about the line 5{x + 1) = — S(x) for both dynamics. 

V. DISCUSSION AND CONCLUSIONS 

In this work we have continued our studies of the microstructure of an unrestricted SOS interface driven far from 
equilibrium by an applied field. Previous studies indicate that different dynamics can lead to important differences 
in the microstructure of the moving interface and that extreme care therefore must be taken in selecting stochastic 
dynamics appropriate for the specific physical system of interest . fAM. Wi ITo| 

For this study we have considered two dynamics that include a local energy barrier representing a transition state 
inserted between the individual Ising or lattice-gas states. Such Arrhenius dynamics, as they are often called, are 
appropriate in kinetic MC simulations of discrete Ising or lattice-gas models in which the discrete states serve as 
approximations for high-probability configurations in an underlying continuous potential. 

The two Arrhenius dynamics that we considered are the commonly used one-step-dynamics lDSD)[24Ll25| and the 
two-step transition dynamics approximation fTDA) . [Til IrH l2t| The OSD belongs to the class of dynamics known as 
soft and the TDA to the class known as hard.pH 

We studied the microstructure and velocity of the SOS interface by kinetic MC simulations and by a non-linear 
mean- field theory developed in previous papers. 0,13 We calculated the interface velocity as a function of the driving 
field, temperature, and angle of the interface relative to the lattice axes. We also studied the local shape of the 
interface in terms of the spin-class populations, the average height of a step, and the probability density for individual 
steps in the interface. The theory predicts significant differences between interfaces moving under hard and soft 
dynamics. For soft dynamics, interface structures should remain independent of the applied field, while there must 
be a clear dependence on the field for hard dynamics. 

For the TDA dynamics we found generally very good agreement between the theoretical predictions and the MC 
simulations. However, for the OSD dynamics we found that, contrary to the theoretical prediction, there is a weak 
but clear dependence of the interface structure on the field. This dependence is manifest in the average stationary 
step height and in the mean stationary class populations, which both show a weak dependence with the field that 
saturates for strong fields. As a consequence of this dependence, the theoretical results for the velocities do not match 
very closely the simulated results for low temperatures, particularly for strong fields and small angles. 

The interfaces moving under the hard TDA dynamics present similar characteristics to those evolving under a hard 
Glauber dynamics. However, we found significant differences between the structure of surfaces evolving under the soft 
OSD dynamics and the soft Glauber dynamics studied in Ref. |{J The velocities under the OSD dynamics present a 
discontinuity at T — that is not observed for the soft Glauber dynamics. More interesting is the existence of strong 
skewness in the OSD model. This indicates that lack of skewness is not a necessary characteristic of soft dynamics, 
as earlier results seemed to suggest. [9l ITo| 

Within the mean-field approximation used here, individual steps of the interface are assumed to be statistically 
independent. Short-range correlations are not taken into account by this approximation. The skewness between the 
spin population on the leading and trailing edges of the interface is a consequence of such short-range correlations. 
For both the TDA and the OSD dynamics, the interfaces undergo a gradual breakdown of up /down symmetry for 
increasing fields, which has also been observed in other examples of driven interfaces. [^J 13^, E3 This breakdown is 
clearly evident for the OSD model, and it is probably the reason why the mean-field approximation misses the weak 
field dependence in the OSD interface structure. 

It is obviously important to note and eventually to understand the discrepancies between the theoretical mean-field 
predictions and the simulation results for the OSD dynamics. However, on a qualitative level the theory predicts 
quite accurately the differences between the structures generated by the two dynamics. The average step height for 
the TDA increases dramatically with increasing field, as accurately predicted by the theory. For the OSD, the step 
height does increase somewhat with H, in contrast to the theoretical prediction of ii-independence, but the increase 
is very small. In comparison with the TDA, the OSD surface remains very smooth, with an average step height well 
below unity. 

As in previous studies, our results indicate strong differences between interfaces moving under different dynamics, 
emphasizing the need for extreme care in selecting the appropriate dynamics for the physical system of interest. Even 
in cases where soft dynamics are the more appropriate choice, as for solidification or adsorption problems where the 
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driving force is a chemical-potential difference. [T2I Il3l 1271 l2gl |29| the results can depend significantly on which soft 
dynamics are chosen. 
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TABLE I: The spin classes in the anisotropic square-lattice SOS model. The first column contains the class labels, jks. The 
second column contains the total field and interaction energy for a spin in each class, E(jks), relative to the energy of the state 
with all spins parallel and H — 0, Eo = — 2(J X + J y ). The third column contains the change in the total system energy resulting 
from reversal of a spin from s to — s, AE(jks). The fourth and fifth columns contain Et — Ei and Ef — Et, respectively. 
The first three classes have nonzero populations in the SOS model, and flipping a spin in any of them preserves the SOS 
configuration. The other two classes (marked f ) also have nonzero populations in the SOS model, but flipping a spin in any 
of them would produce an overhang or a bubble and is therefore forbidden. Not that in this table, s represents the spin value 
before the spin flip. 



Class, jks 


E(jks) - Eo 


AE(jks) 


Et — Ei 


Ef — Et 


01s 


-sH + 2J y 


2sH + 4 J x 


sH + 2J X + U 


sH + 2J X - U 


lis 


-sH + 2{J X + Jy) 


2sH 


sH + U 


sH -U 


21s 


-sH + 2{2J X + J y ) 


2sH - AJ X 


sH - 2J X + U 


sH - 2J X - U 


10s f 


-sH + 2J x 


2sH + 4 J H 


sH + 2Jy + U 


sH + 2J y -U 


20s f 


-sH + 4J x 


2sH - 4(J X - J y ) 


sH - 2{J X -J y ) + U 


sH - 2(J X -J y )-U 



TABLE II: The mean populations for the spin classes of the SOS interface, with the corresponding contributions to the 
interface velocity under the TDA and the OSD dynamics. The first column contains the class labels, jks. The second column 
contains the mean spin-class populations for general tilt angle <f), with cosh7(<£) from Eq. J2j. The third and fourth columns 
contain the contributions to the mean interface velocity in the y direction from spins in classes jk— and jk+, Eq. 10, using the 
SOS-preserving TDA and OSD dynamics respectively. For the TDA dynamics, X — X(T, H) is given by Eq. If 21 . and for the 
OSD dynamics X is independent of H, and is given by Eq. 1131 . In the third column A = [cosh (2/3 J x ) cosh (/3H) + cosh f3U] , 
and B = sinh (2/3 J x ) sinh (J3H). 



Class, jks 


(n(jks)) 


<«y(jfc)>TDA 


(Vy(jk)) OSD 


01s 


1-2X cosh~,(4>) + X'' ! 


-2/3 J x Asinh(pH)+B cosh((3/J) 




(1-X a ) 2 




lis 


2X\(1+X' Z ) cosh j(4>)-2X] 


sinh(/3/J) 


e - /3l/ 2sinh (J3H) 


(l-X")" 


cosh(,37r/) + cosh f3U 


21s 


X' ! [1-2X cosl^OJ+X^] 


2/3 J x Asinh(/3/J)-flcosh(^H) 


e -W-2J*) 2sinh {(jH) 
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FIG. 1: A short segment of an SOS interface y = h(x) between a positively magnetized phase (or "solid" phase in the lattice-gas 
picture) below and a negative (or "fluid") phase above. The step heights are 5(x) = h(x + 1/2) — h(x — 1/2). Interface sites 
representative of the different SOS spin classes (see Tableland Table |HJ are marked with the notation jks explained in the 
text. Sites in the uniform bulk phases are 00— and 00+. This interface was randomly generated with a symmetric step-height 
distribution, corresponding to (f> = 0. From Ref. 0. 



(Ej + E f ) 
2 




initial transition final 
state state state 



FIG. 2: Schematic picture of the transition barrier in the symmetric Butler- Volmer approximation, used to calculate the TDA 
and OSD transition rates. After Ref. I23I 
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FIG. 3: (Color online) MC (data points) and analytical (solid lines) results for the stationary single-step pdf calculated with 
the TDA dynamics, shown on a logarithmic scale vs 8, for the values of H/ J given in the legend, (a) T = 0.2T C . (b) T — 0.6T C . 
The symbols (and colors) have the same interpretations in (a) and (b). 




FIG. 4: (Color online) The stationary single-step pdf calculated with the OSD dynamics, shown on a logarithmic scale vs S. 
The data points indicate MC results and the straight lines are the theoretical predicted values (independent of H). The values 
of H/J are given in the legend, (a) T = 0.2T C . (b) T = 0.6T C . The symbols and colors have the same interpretations in (a) 
and (b). Note the very different scales from Fig.|H] 
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FIG. 5: Average stationary step height (\S\) vs H for <f>=0 at T=0.2T C and 0.6T C . The curves represent the theoretical results. 
The MC data were obtained directly by summation over the simulated single-step pdfs (asterisks and crosses) and from the 
probability of zero step height (circles and squares). See the text for details. Curve with circles and asterisks: T = 0.2T C . Curve 
with squares and crosses : T = 0.6T C . (a) TDA dynamics, shown on a logarithmic vertical scale, (b) OSD dynamics, shown on 
a linear vertical scale. In this and all the following figures, the statistical uncertainty is much smaller than the symbol size. 




FIG. 6: The average stationary normal interface velocity (v±) vs H for <fi = 0. The MC results are shown as data points, 
circles for T = 0.2T C and squares for T = 0.6T C , and the theoretical results as solid curves, (a) TDA, shown on a linear vertical 
scale, (b) OSD, shown on a logarithmic vertical scale. 
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FIG. 7: The average stationary normal interface velocity (v±) vs tan</>, calculated with the TDA dynamics, for H/J — 0.1 
(up triangles), 0.5 (circles), 1 (squares), 1.5 (diamonds), and 2.0 (left triangles). The symbols represent MC data, and the solid 
curves analytical results, (a) T = 0.2T C . (b) T = 0.6T C . 




tan((p) 



<] 



< > O O O O ^gjxc 



< 



< 




□ □□□□□ e-^^^^^^^^ 



OOOOOOOOOOO Q-O-Q-Q-gH^e-Q-O-Q 



Aaaa^aaaaaaaaa^aaaaaz \ 

u n o.2 n zl n a n s i 



0.4 0.6 
tan((p) 



0.8 



FIG. 8: The average stationary normal interface velocity (v±) vs tan</>, calculated with the OSD dynamics, for H/J — 0.1 (up 
triangles), 0.5 (circles), 1 (squares), 1.5 (diamonds), and 2.0 (left triangles) . The symbols represent MC data, and the solid 
curves analytical results, (a) T = 0.2T C . (b) T — 0.6T C . Note the difference in the vertical scales. 
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T/Tc T/Tc 

FIG. 9: (Color online) The average stationary normal interface velocity (v±) vs T for = and H/J between 0.5 and 3.5. 
MC data are represented by data points and analytical results by solid curves. From below to above, the values of H/J are 0.5 
(circles), 1.0 (squares), 1.5 (diamonds), 2.0 (up triangles), 2.5 (left triangles), 3.0 (down triangles), and 3.5 (right triangles). 
Online, the colors of the curves and symbols match, (a) TDA, on a linear vertical scale, (b) OSD, on a logarithmic vertical 
scale. 




FIG. 10: Mean stationary class populations (n(jks)) vs H/J for <j> — 0, calculated for the TDA dynamics. The simulation 
results are indicated by symbols, and the analytic approximations by solid curves, (a) T — 0.2T C . (b) T = 0.6T C . The symbols 
have the same interpretations in (a) and (b). Note the different vertical scales in the two parts. 
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FIG. 11: Mean stationary class populations (n(jks)) vs H/J for = calculated for the OSD dynamics. The simulation 
results are indicated by the symbols, and the straight lines indicate the theoretical predicted values (independent of H). (a) 
T — 0.2T C . (b) T = 0.6T C . The symbols have the same interpretations in (a) and (b). Note the different vertical scales in the 
two parts. 
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FIG. 12: The two relative skewness parameters p (circles) and e (squares), defined in Eqs. 1 1H and IJ5L respectively. The 
parameters are shown vs H for <f> = 0, at T = 0.2T C (empty symbols) and at T = 0.6T C (filled symbols), (a) TDA. (b) OSD. 
Note the different vertical scales in the two parts. 
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FIG. 14: Contour plots of log 10 p [S(x), S(x + 1)] for = at T = 0.6T C for the OSD dynamics, (a) H/J = 0. (b) H/J = 1.0. 
(c) J//J = 2.0. (d) H/J = 2.5. See discussion in the text. 



18 



